clear; clc; close all;

k = 3;
c1 = 1e-9; c2 = 0; c3 = 0;

eta = 0;  % 固定 eta

h = 0.01;
T = 600;
transient =550;
N = T / h;
N_trans = transient / h;

mu_list = linspace(-3*pi, 3*pi, 300);

Y_peaks = cell(length(mu_list), 1);

parfor i = 1:length(mu_list)
    mu = mu_list(i);
    
    x = [0; 0; 0];
    Y_data = zeros(N,1);
    
    for t = 1:N
        x = RK4_step(@(x) reduced_system_ode(x,k,c1,c2,c3,mu,eta), x, h);
        Y_data(t) = x(2);
    end
    
    Y_steady = Y_data(N_trans+1:end);
    peaks = findpeaks(Y_steady);
    Y_peaks{i} = peaks;
end

figure; hold on; box on;
for i = 1:length(mu_list)
    mu = mu_list(i);
    yvals = Y_peaks{i};
    xvals = mu * ones(size(yvals));
    plot(xvals, yvals, 'k.', 'MarkerSize', 2);
end

xlabel('\mu');
ylabel('Y_{max}');
xticks(-3*pi : pi : 3*pi);
xticklabels({'-3\pi','-2\pi','-\pi','0','\pi','2\pi','3\pi'});
grid on;

%% 下面函数保持不变

function dx = reduced_system_ode(x, k, c1, c2, c3, mu, eta)
    X = x(1);
    Y = x(2);
    Z = x(3);
    dx = zeros(3,1);
    
    dx(1) = Y + Z - k * sin(Y + eta) + k * sin(eta) + c1;
    dx(2) = -X + Z + c2;
    dx(3) = -X - Z + k * sin(X + mu) - k * sin(mu) + c3;
end

function x_next = RK4_step(f, x, h)
    k1 = f(x);
    k2 = f(x + h/2 * k1);
    k3 = f(x + h/2 * k2);
    k4 = f(x + h * k3);
    x_next = x + h/6*(k1 + 2*k2 + 2*k3 + k4);
end
